##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(ggplot2)
library(ggsci)
library(ggpubr)
library(patchwork)
library("scales")
library(parallel)

##########################################################################################

option_list <- list(
    make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--sample_list_public_file"), type = "character"),
    make_option(c("--rsem_file"), type = "character"),
    make_option(c("--maf_public_file"), type = "character"),
    make_option(c("--gene"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"

    sample_list_file <- "~/20220915_gastric_multiple/dna_combinePublic/config/tumor_normal.class.list"
    sample_list_public_file <- paste(work_dir,"/public_ref/combine/MutationInfo.combine.tsv",sep="")

    rsem_file <- "~/20220915_gastric_multiple/dna_combinePublic/mRNA/CombineTMM.DNAUse.NJMU_TCGA.tsv"
    maf_public_file <- paste(work_dir,"/maf_public/All_use.addVAF.maf",sep="")
    
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/images/expression"
    gene <- "MUC6"

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
out_path <- opt$out_path
rsem_file <- opt$rsem_file
sample_list_public_file <- opt$sample_list_public_file
maf_public_file <- opt$maf_public_file
gene <- opt$gene

dir.create( out_path , recursive = T )

##########################################################################################

info <- data.frame(fread(sample_list_file))
info_public <- data.frame(fread(sample_list_public_file))
dat_expression <- data.frame(fread(rsem_file))
dat_maf_public <- data.frame(fread( maf_public_file ))
use_gene <- gene

##########################################################################################

Variant_Types <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")
dat_maf_public <- subset( dat_maf_public , Hugo_Symbol == gene & Variant_Classification %in% Variant_Types )
colnames(dat_expression) <- gsub( "[.]" , "-" , colnames(dat_expression) )

##########################################################################################

info_public <- subset( info_public , From != "NJMU" )
info_public$ID <- info_public$Tumor
info_public$Class_sub <- info_public$Class
info_use <- rbind( info_public[,c( "ID" , "Tumor" , "Class" , "Class_sub" , "From")] , info[,c( "ID" , "Tumor" , "Class" , "Class_sub" , "From")] )

mutTumor <- unique(dat_maf_public$Tumor_Sample_Barcode)
info_mut <- subset( info_use , Tumor %in% mutTumor )
info_mut <- paste0(info_mut$ID , "_" , info_mut$Class_sub)
info_mut <- info_mut[info_mut %in% colnames(dat_expression)]

info_wild <- subset( info_use , !(Tumor %in% mutTumor) )
info_wild <- paste0(info_wild$ID , "_" , info_wild$Class_sub)
info_wild <- info_wild[info_wild %in% colnames(dat_expression)]

##########################################################################################

getTMM <- function(exp_use = exp_use , type = type){
    sample <- sapply( strsplit(names(exp_use) , "_") , "[" , 1)
    class <- sapply( strsplit(names(exp_use) , "_") , "[" , 2)
    class <- sapply( strsplit(class , "-") , "[" , 1)
    tmm <- as.numeric(exp_use)
    tmp_dat <- data.frame( Sample = sample , Class = class , TMM = tmm )
    ## 合并一个人的同一样本
    tmp_dat <- tmp_dat %>%
    group_by( Sample , Class ) %>%
    summarize( TMM = median(TMM) )
    tmp_dat <- data.frame(tmp_dat)
    tmp_dat$Type <- type

    return(tmp_dat)
}

##########################################################################################
## 只看肠化里面

result <- Reduce(function(x,y)bind_rows(x,y) , mclapply( unique(dat_expression$gene_id) , function(gene){
    print(gene)

    tmp_exp <- subset( dat_expression , gene_id == gene )
    
    ## 突变型样本的表达
    type <- paste0(use_gene,"_Mut")
    use_sample <- info_mut
    exp_use <- tmp_exp[use_sample]
    tmm_mut <- getTMM(exp_use = exp_use , type = type)

    ## 野生型样本的表达
    type <- paste0(use_gene,"_Wild")
    use_sample <- info_wild
    exp_use <- tmp_exp[use_sample]
    tmm_wild <- getTMM(exp_use = exp_use , type = type)

    ## 合并
    tmm_combine <- rbind( tmm_wild , tmm_mut )
    tmm_combine$gene <- gene

    ## 只提取IM
    tmm_combine <- subset( tmm_combine , Class == "IM" )

    mut <- subset( tmm_combine , Type == paste0(use_gene,"_Mut") )
    wild <- subset( tmm_combine , Type == paste0(use_gene,"_Wild") )
    p_value <- wilcox.test(mut$TMM , wild$TMM)$p.value
    tmp_diff <- data.frame( gene = gene , median_tpm_mut = median(mut$TMM) , median_tpm_wild = median(wild$TMM) , p = p_value )
    return( tmp_diff )

},mc.cores=30))


result$FoldChange <- result$median_tpm_mut/result$median_tpm_wild
result$fdr <- p.adjust( result$p , method = "fdr" )

out_name <- paste0( out_path , "/MUC6_AllDiff.MutVsWild.IM.tsv" )
write.table(result , out_name , row.names = F , sep = "\t" , quote = F)

